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ABSTRACT 

Spatial dependence in the statistics of redshifted 21cm fluctuations promises to 
provide the most powerful probe of the reionisation epoch. In this paper we consider 
the second and third moments of the redshifted 21cm intensity distribution using a 
simple model that accounts for galaxy bias during the reionisation process. We demon- 
strate that skewness in redshifted 21cm maps should be substantial throughout the 
reionisation epoch and on all angular scales, owing to the effects of galaxy bias which 
leads to early reionisation in over-dense regions of the IGM. The variance (or power- 
spectrum) of 21cm fluctuations will exhibit a minimum in redshift part way through 
the reionisation process, when the global ionisation fraction is around 50%. This min- 
imum is generic, and is due to the transition from 21cm intensity being dominated 
by over-dense too under-dense regions as reionisation progresses. We show that the 
details of the reionisation history, including the presence of radiative feedback are 
encoded in the evolution of the auto-correlation and skewness functions with red- 
shift and mean IGM neutral fraction. The amplitudes of fluctuations are particularly 
sensitive to the masses of ionising sources, and vary by an order of magnitude for 
astrophysically plausible models. We discuss the detection of skewness by first gen- 
eration instruments, and conclude that the Mileura Widefield Array-Low Frequency 
Demonstrator will have sufficient sensitivity to detect skewness on a range of angular 
scales at redshifts near the end of reionisation, while a subsequent instrument of 10 
times the collecting area could map out the evolution of skewness in detail. The ob- 
servation of a minimum in variance during the reionisation history, and the detection 
of skewness would both provide important confirmation of the cosmological origin of 
redshifted 21cm intensity fluctuations. 

Key words: cosmology: diffuse radiation, large scale structure, theory - galaxies: 
high redshift, inter-galactic medium 



1 INTRODUCTION 

The primordial density field as revealed in maps of the Cos- 
mic Microwave Background may be described as a Gaus- 
sian random field. In other words, the distribution of over- 
densities in the Fourier decomposition of the field is Gaus- 
sian at all spatial frequencies, so that the statistics of the 
field may be described by a single number (the variance) as 
a function of spatial scale. This scale dependent variance is 
referred to as the power-spectrum. On large enough scales 
we might also expect the statistics for 21cm emission to be 
Gaussian, since over-dense regions have more hydrogen, and 
therefore a larger 21cm intensity. On small scales this Gaus- 
sianity will be broken by the presence of HII regions. Thus 



much attention has focused on the prospects of measuring 
the power-spectrum of 21cm emission, as well as the use of 
the power-spectrum to probe the astrophysics of reionisa- 
tion (e.g. Zaldariaga et al. 2004; Morales & Hewitt 2004; 
Furlanetto et al. 2004). 

In regions that are over-dense, galaxies will be over- 
abundant for two reasons: first because there is more mate- 
rial per unit volume to make galaxies, and second because 
small-scale fluctuations need to be of lower amplitude to 
form a galaxy when embedded in a larger-scale over-density 
(the so-called galaxy bias; see Mo & White 1996). Regard- 
ing reionisation of the intergalactic medium (IGM) , the first 
effect will result in a larger density of ionising sources. How- 
ever this larger density will be compensated by the increased 
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density of gas to be ionised. Furthermore, the increase in the 
recombination rate in over-dense regions will be counter- 
acted by the galaxy bias in over-dense regions. The process 
of reionisation also contains several layers of feedback. Ra- 
diative feedback heats the IGM and results in the suppres- 
sion of low-mass galaxy formation (Efstathiou, 1992; Thoul 
& Weinberg 1996; Quinn et al. 1996; Dijkstra et al. 2004). 
This delays the completion of reionisation by lowering the 
local star formation rate, but the effect is counteracted in 
over-dense regions by the biased formation of massive galax- 
ies. The radiation feedback may therefore be more important 
in low-density regions where small galaxies contribute more 
significantly to the ionising flux. 

Wyithe & Loeb (2007) have modeled the density depen- 
dent reionisation process using a semi-analytic model that 
incorporates the features described above, and so captures 
the important physical processes. This model demonstrates 
that galaxy bias leads to enhanced reionisation in over-dense 
regions. Because this bias operates on the exponential tail of 
the Press-Schechter (1974) mass function, the enhancement 
of ionisation is not linear with over-density. As a result, 
galaxy bias leads to a non-Gaussian distribution of 21cm 
brightness temperature in intensity maps of reionisation. 
This effect was discussed by Lidz et al. (2006) using detailed 
numerical simulations. Lidz et al. (2006) found that the con- 
tribution of a third order term in the power-spectrum will 
be comparable to the second order term at an epoch where 
the IGM is approximately 50% ionised. While comparison 
with data may require detailed numerical simulations of the 
sort described in Lidz et al. (2006), computationally cheap 
semi-analytic models are useful for exploring which physi- 
cal parameters will be probed by the detection of angular 
fluctuations in the 21cm intensity. In this paper we explore 
the dependence of observables on the values of parameters 
of interest like the clumping of the IGM and the minimum 
mass for galaxy formation. We also discuss the prospects 
for detection of the third moment (skewness) with first and 
second generation low-frequency arrays. 

In § [5] and § [3] we describe our model for density de- 
pendent reionisation, and the resulting redshifted 21cm in- 
tensity probability distribution (including calculation of its 
second and third moments). We describe the features of 
the auto-correlations and skewness functions that corre- 
spond to different physical parameters that are important 
for the reionisation epoch. We then discuss the measure- 
ment of these moments using first and second generation 
low-frequency arrays (in particular the Mileura Widefield 
Array-Low Frequency Demonstrator) in §0] before conclud- 
ing in § [5] Throughout the paper we adopt the set of cos- 
mological parameters determined by WMAP (Spergel et al. 
2006) for a flat ACDM universe. 



2 DENSITY DEPENDENT MODEL OF 
REIONISATION 

In this paper we compute the relation between the local 
over-density and the brightness temperature of redshifted 
21cm emission using the model described in Wyithe & 
Loeb (2007). Here we summarise the main features of the 
model and refer the reader to that paper for more details. 
The evolution of the ionisation fraction by mass Qs,r 



of a particular region of scale R with over-density S (at ob- 
served redshift z bs) may be written as 



dQ s 



dt 



0.76 



dF col (S,R,z,M io 



(l-Qs,n 



dt 

dF col (5,J?,2,M min ) 



dt 



(1) 



where ]Vi on is the number of photons entering the IGM 
per baryon in galaxies, Ob is the case-B recombination co- 
efficient, C is the clumping factor (which we assume, for 
simplicity, to be constant), and D(z) is the growth factor 
between redshift z and the present time. The production 
rate of ionising photons in neutral regions is assumed to be 
proportional to the collapsed fraction F co \ of mass in halos 
above the minimum thresholds in neutral (M m t n ), and in 
ionised (Mi on ) regions. For our fiducial model, we assume 
Afmin to correspond to a virial temperature of 10 4 K, repre- 
senting the hydrogen cooling threshold, and Mi on to corre- 
spond to a virial temperature of 10 K, representing the mass 
below which infall is suppressed from an ionised IGM (Dijk- 
stra et al. 2004). In a region of co-moving radius R and mean 
over-density 5(z) = 8D(z) / ' D{z ^ s ) [specified at redshift z 
instead of the usual 2 = 0], the relevant collapsed fraction 
is obtained from the extended Press-Schechter (1974) model 
(Bond et al. 1991) as 



F co i(S, R, z) = erfc 



S(z) 



O-galJ 



HR)Y 



(2) 



where erfc(a;) is the error function, a(R) is the variance of 
the density field smoothed on a scale R, and cr ga i is the 
variance of the density field smoothed on a scale R sa i , corre- 
sponding to a mass scale of M m i n or Mi on (both evaluated at 
redshift z rather than at z = 0). In this expression, the criti- 
cal linear over-density for the collapse of a spherical top-hat 
density perturbation is 8 C ~ 1.69. 

The model assumes that on large (linear-regime) scales 
most ionising photons are absorbed locally, so that the ion- 
isation of a region is caused by nearby ionisation sources. 
This assumption is certainly justified during the early stages 
of reionisation, when the mean free path for ionising pho- 
tons is short. However even later in the reionisation process, 
the mean free-path always remains smaller than the charac- 
teristic HII bubble size (it could be smaller if mini-halos or 
pockets of residual HI block ionising photons between the 
sources and the edge of the HII region.) Our local ionisation 
assumption is therefore valid as long as the characteristic 
bubble size is smaller than the spatial scale of the correla- 
tions we consider. In this paper we consider scales of 1-25 
arc-minutes which will be of interest to upcoming 21cm ex- 
periments, corresponding to co-moving scales of 2-50 Mpc. 
These scales are larger than the mean-free-path of ionising 
photons near the end of reionisation (Fan et al. 2006). 

As an example, we find the value of iVi on that yields 
overlap of ionised regions at the mean density IGM by z ~ 6 
(White et al. 2003). Equation |T} may be integrated as a 
function of 8. At a specified redshift, this yields the filling 
fraction of ionised regions within the IGM on various scales 
R as a function of over-density. We may then calculate the 
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Figure 1. The difference in brightness temperature between over 
or under-dense regions, and the IGM at mean density as a func- 
tion of redshift. The lines represent regions at redshift z that are 
—3cr, —la, —let, Oct, icr, 2<r or 3cr fluctuations in the smoothed 
density field. We assumed a clumping factor of C = 10 and an 
angular scale of 10'. 



corresponding 21cm brightness temperature contrast 



T{5, R) = 22mK(l - Q S , R ) \1 + -5 J , (3) 

where the pre- factor of 4/3 on the over-density refers to 
the spherically averaged enhancement of the brightness tem- 
perature due to peculiar velocities in over-dense regions 
(Bharadwaj & Ali 2005; Barkana & Loeb 2005). 

The predicted evolution of the brightness temperature 
of regions of IGM at fixed angular scale is plotted in Figure[T] 
For clarity we have plotted T — T(S = 0), i.e. the difference 
in brightness temperature between the over or under-dense 
region, and the IGM at mean density. Each line in the fig- 
ure represents a region that is a — 3a, —2a, —la, 0<r, la, 
2a or 3ct fluctuation in the smoothed density field at red- 
shift z. Here we assumed a clumping factor of C = 10 and 
an angular scale of 10', where 9 — R/d\(z) and d\ is the 
angular diameter distance. At high redshift the curves are 
symmetric about zero, as expected for linear fluctuations 
in the primordial density field. As reionisation proceeds the 
brightness temperature (relative to the mean IGM) is ini- 
tially quite insensitive to redshift. However at later times 
the brightness temperature of over-dense regions decreases, 
while the brightness temperature of under-dense regions in- 
creases (again when compared to the mean IGM). Finally, 
reionisation is completed in over-dense regions of the IGM 
first, as may be seen in the 3a fluctuation which is reionised 
at z ~ 7, prior to the mean IGM at z = 6. Therefore at 
some point during the reionisation history there must be a 
transition between bright patches of 21cm emission being 
observed in over-dense regions, too bright patches being ob- 
served in under-dense regions. This transition results in an 
epoch of nearly uniform brightness temperature (uniform 
when smoothed on an angular scale larger than the typical 
HII region) of the IGM. Figure [1] shows that in our models 
this transition occurs near a redshift of z ~ 8 — 9. The epoch 
of near uniform temperature across the IGM has interesting 
consequences for the variance in smoothed 21cm intensity 
maps. This topic will be discussed in the following section. 

Figure [1] also shows that reionisation skews the bright- 



ness temperature towards low values, owing to the dominant 
effect of galaxy bias. A point of interest here is that once 
reionisation begins, this results in the brightness tempera- 
ture of both over- and under-dense regions being lower than 
the brightness temperature of the IGM at mean density (we 
note that the mean is not found to depend on angular scale 
as required for consistency). Since the underlying probabil- 
ity distribution of over-densities dPj dS is known (a Gaussian 
of variance aD), we may compute the observed probability 
distribution for T, 



dP dP 



85_ 
dT 



(4) 



This distribution corresponds to fluctuations in redshifted 
21cm emission on angular scales 9. Example distributions 
are shown in Figures [2] and [3] 

Firstly we discuss the evolution of the distribution with 
redshift at fixed angular scale. In Figure [21 distributions are 
shown assuming a clumping factor of C = 10 and an an- 
gular scale of 10', at redshifts of z — 6.5, z — 7, z = 8.5, 
2 = 10, z = 13 and z — 22. The large symbols in this 
figure show the value of T — T(8 = 0) corresponding to re- 
gions that are —2a, —la, 0a, la and 2a fluctuations in the 
smoothed density field respectively. The lower right hand 
panel clearly shows the Gaussian distribution at high red- 
shift, with evenly spaced values of T corresponding to the 
±lcr and ±2cr points. At slightly lower redshift, the lower 
central panel shows that the Gaussian has become a little 
wider due to linear growth, with the high T tail slightly 
truncated (due to the onset of skewness), as may be seen 
by comparing the relative likelihoods for temperature corre- 
sponding to the positive and negative over-density fluctua- 
tions. At a redshift of ~ 10 in this example (lower left panel), 
one gets a narrower and very skewed distribution. At this 
point the positive over-densities still correspond to larger 
than average brightness temperature, but the difference in 
brightness temperature between regions of mean density and 
+2(7 fluctuations is small. The upper right panel shows the 
distribution at z — 8.5 (the minimum). At this redshift, 
the distribution is both very asymmetric and narrow, but in 
addition, the monotonic dependence of brightness temper- 
ature with over-density is broken. Indeed, the mean IGM 
density corresponds to the largest brightness temperatures, 
while the brightness temperatures of +la and —la, and 
of +2(7 and — 2cr density fluctuations are almost the same. 
At lower redshifts (upper-central and upper-left panels) the 
width and skewness of the distributions both increase, and 
a monotonic relation between the over-density and bright- 
ness temperature is restored, though the order is reversed 
relative to high redshift. 

In Figure [3] distributions are shown assuming a clump- 
ing factor of C = 10 and a redshift of z = 7, at angular- 
scales of 5', 10', 15' and 20' and z = 7. The distribution is 
valid on angular scales sufficiently large that no region on 
that scale is fully ionised. On small enough scales the distri- 
bution must be described using a model that calculates the 
size distribution of individual HII regions (e.g. Furlanetto et 
al. 2004). However as may be seen from Figure [3] our model 
describes the distribution dP/dT for 9 > 5' at z = 7, and 
for all angles of interest to up-coming 21cm experiments at 
higher redshifts. At large angular scales (lower-right panel) 
the distribution is narrow and close to Gaussian. By compar- 
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T-T(S=0) T-T(8=0) T-T(8=0) 



Figure 2. The probability distribution of brightness temperature dP/dT. The horizontal axis shows the difference in brightness temper- 
ature between over or under-dense regions, and the IGM at mean density [T — T(<5 = 0)]. Distributions are shown assuming a clumping 
factor of C = 10 and an angular scale of 10', at redshifts of z = 6.5, z = 7, z = 8.5, z = 10, z = 13 and z = 22. The large symbols 
show the value of T — T{6 = 0) corresponding to regions that are — 2a (filled circle), — la (filled square), Oct (filled diamond), la (empty 
square) and 2a (empty circle) fluctuations in the density field and correspond to lines in Figure [T] 



ing the brightness temperatures corresponding to the Tier 
and ±2<r regions of IGM with the redshift dependent distri- 
butions in Figure [5] we see that even though the distribu- 
tion is still close to Gaussian it already includes the effects 
of galaxy bias, since the brighter than average regions are 
formed by under-dense regions of IGM. As the angular scale 
is decreased the width of the distributions increases, and the 
skewness towards small T increases. 

Figures [2] and [3] explicitly demonstrate the introduction 
of skewness into the distribution of T(9). The skewness of 
21cm intensity maps, and its dependence on model parame- 
ters will be discussed in more detail in the following section. 



3 THE VARIANCE AND SKEWNESS OF 
21CM MAPS 

In this section we calculate the second and third moments of 
the probability distributions of T in redshifted 21cm inten- 
sity maps, as functions of angular scale and redshift. Given 
the structure of our semi-analytic model for reionisation, it 
is natural to discuss the moments of the real space intensity 
fluctuations smoothed within top-hat window functions of 
angular radius 9. 

The second moment of these distributions corresponds 
to the auto-correlation function of brightness temperature 
smoothed on an angular radius 9: 



<(T-<T» 2 ) = 
1 



2na{R) 



d8(T(8,R)-{T))e 



where 



(5) 



(6) 



and a{R) is the variance of the density field (at redshift z) 
smoothed on a scale R. 

A measure of the departure of the statistics from those 
of a Gaussian random field is provided by the third moment 
of the distribution dP/ dT. We refer to the third moment as 
the skewness function of brightness temperature smoothed 
with top-hat windows of angular radius 9. 



{(T-{T}Y) = 
1 



2na(R) 



! 



d8(T(8,R)-{T)y e 



(7) 



Examples of auto-correlation and skewness functions 
are shown in the left and right hand panels of Figure [4] as- 
suming clumping factors of C = 10, C — 20 and (7 = 2. 
The auto-correlation function is plotted as a function of 9 
(at z = 6.5, upper panel), as a function of mean IGM neu- 
tral fraction xhi (at 9 — 10', central panel), and as a func- 
tion of redshift (at 9 — 10', lower panel). In the left and 
right hand panels of Figure [5] we show examples of auto- 
correlation and skewness functions for four commonly dis- 
cussed scenarios of galaxy formation in a neutral and ionised 
IGM. The four cases presented in Figure [5] show results for 
the following combinations of (T m in, lion): (10 4 K, 10 5 K) cor- 
responding to our fiducial model with star-formation down 
to the hydrogen cooling threshold in neutral regions plus 
radiative feedback; (10 4 K, 10 4 K) corresponding to a model 
with no radiative feedback; and (10 5 K, 10 5 K) which could 
correspond to a model where feedback within the galaxy 
(e.g. from supernovae Dekel & Woo 2003) limits the star 
formation efficiency in low mass galaxies; and (250K, 10 5 K) 
corresponding to a model where molecular hydrogen cooling 
allows star formation to occur in mini-halos below the hy- 
drogen cooling threshold. The clumping factor was assumed 
to be C = 10 in each case. In all examples presented in Fig- 
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Figure 3. The probability distribution of brightness temperature dP/dT. The horizontal axis shows the difference in brightness temper- 
ature between over or under-dense regions, and the IGM at mean density [T — T(& = 0)]. Distributions are shown assuming a clumping 
factor of C = 10 and a redshift of z = 7, at angular scales of 5', 10', 15' and 20'. The large symbols show the value of T — T(<5 = 0) 
corresponding to regions that are —2a, — la, Oct, Ict and 2ct fluctuations in the density field respectively. 



ures l4!5l we have set the efficiency of the ionising sources so 
as to achieve reionisation in the mean IGM at z — 6. 

We find that the amplitude of the auto-correlation func- 
tion decreases with angular scale as expected. The ampli- 
tude of the skewness also decreases with angular scale. How- 
ever interestingly, the amplitudes of the auto-correlation and 
skewness functions do not vary monotonically with redshift 
or average neutral fraction. In particular, both the variance 
and skewness of the 21cm intensity field have minima dur- 
ing the reionisation history. At high redshift the skewness 
goes to zero, since the 21cm map prior to galaxy formation 
must follow a Gaussian distribution corresponding to the 
primordial density fluctuations. 

This unintuitive behavior was alluded to in Figure [T] 
and can be understood as follows. Early in the reionisation 
history, before galaxies have had much influence on the ioni- 
sation state of the IGM, the variance is governed by the pri- 
mordial power-spectrum of fluctuations in the density field, 
and the skewness is zero. During this early time, the over- 
dense regions are brighter due to the higher density of neu- 
tral IGM there. As galaxies modify the ionisation state of 
the IGM, the over-dense regions get ionised first, breaking 
the Gaussianity of the 21cm map. Therefore, once reioni- 
sation is underway, the over-dense regions of IGM become 
fainter in 21cm emission relative to under-dense regions. As 
was shown in Figure [1] this results in a transition between 
bright patches of 21cm emission being observed in over- 
dense regions, too bright patches being observed in under- 
dense regions. This transition yields an epoch of nearly uni- 
form brightness temperature of the IGM when smoothed at 
a certain scale, and hence minima in the auto-correlation 
and skewness functions. At later times, both the variance 
and the skewness increase as reionisation is completed. 

To illustrate this complex behavior in more detail, we 
plot contours of the auto-correlation and skewness in Fig- 
ure [6] as functions of redshift and angular scale (upper and 



central panels). For comparison we also plot the evolution 
of the average neutral fraction with redshift (lower panels). 
Two cases are shown, corresponding to the two models in 
Figure [5] that include radiative feedback, i.e. (T m i n , Tion) = 
(10 4 K,10 5 K) [left panels] and (T min ,T ion ) = (250K, 10 5 K) 
[right panels]. Comparison of the upper and central pan- 
els with the lower panels in Figure [6] shows that minima 
in both the auto-correlation and skewness occur when the 
ionised fraction is around Q ~ 0.5 — 0.7, near a redshift of 
z ~ 8 — 9. The minimum occurs at higher redshifts for fluc- 
tuations measured at smaller angular scales. Comparison of 
the central and lower panels of Figure [6] show that skewness 
appears in the 21cm map at the redshift where the reioni- 
sation of the IGM begins (at z ~ 20 in this model). How- 
ever the upper panels show non-zero fluctuations in 21cm 
emission at earlier times, which are driven by the density 
fluctuations in the IGM. Note that the two models yield 
identical fluctuation statistics at z ~ 20. Finally, Figure [6] 
shows that the skewness becomes larger in comparison to 
the auto-correlation amplitude as reionisation proceeds. 

In Figures [4] and [5] we demonstrated the dependence of 
the auto-correlation and skewness functions on model pa- 
rameters as a function of angle or redshift. Figure [6] demon- 
strates this dependence more-fully. While the general trends 
and overall shape are similar in both cases, there are sub- 
stantial quantitative differences. These plots demonstrate 
the rich detail on the reionisation history that can be found 
by considering both the auto-correlation and skewness func- 
tions and their evolution. In particular, the amplitudes of 
the auto-correlation and skewness functions, as well as the 
epoch and depth of their minima in redshift are sensitive 
to the details of the reionisation history. The shape of the 
angular power-spectrum has been shown to encode informa- 
tion on the topology of the reionisation process. Our results 
show that the evolution of the fluctuation statistics in 21cm 
intensity maps will provide complementary information on 
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the astrophysics of the reionisation epoch. This is the topic 
of the next section. 

3.1 Variation of auto-correlation and skewness 
with model parameters 

The free parameters that govern the details of the reioni- 
sation history described in equation JT} are the clumping 
factor C and the minimum masses for galaxy formation in a 
neutral (Mmin) and reionised IGM (Mi on ). These parameters 
encapsulate the astrophysics of the IGM, and of the connec- 
tion between the IGM and ionising sources respectively. In 
this section we discuss the variation of the auto-correlation 
and skewness functions with these parameters. 

In Figure [4] we demonstrated the variation of the angu- 
lar auto-correlation and skewness functions with clumping 
of the IGM. A larger clumping factor results in an increased 
recombination rate in over-dense regions. At low redshift, 
where over-dense regions have lower than average 21cm in- 



tensity, this increased recombination rate reduces the effect 
of galaxy bias on the variation of 21cm intensity with over- 
density by delaying the reionisation of over-dense regions. 
Thus the observational effect of a large clumping factor in 
the IGM is to reduce both the variance and the skewness of 
the distribution of 21cm intensity fluctuations. Conversely at 
higher redshifts (where over-dense regions have higher than 
average 21cm intensity), the increased recombination rate 
results in an increased amplitude for the auto-correlation 
and skewness functions. 

In Figure [S] we demonstrate the dependence of the angu- 
lar auto-correlation and skewness functions on the minimum 
mass for galaxy formation in a neutral and ionised IGM. The 
four cases shown correspond to commonly discussed possi- 
bilities for the star formation history at high redshift. Our 
fiducial model includes radiation feedback by increasing the 
virial temperature corresponding to the minimum mass from 
?mm to Hon following reionisation of a region. As a result, 
once part of a region becomes ionised, the photon produc- 
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Figure 5. Left: Examples of auto-correlation functions, assuming sets of (T min ,T ion ) = (10 4 K, 10 5 K), (10 4 K, 10 4 K), (10 5 K, 10 5 K) and 
(250K, 10 5 K). The clumping factor was chosen to be C = 10. The auto-correlation is plotted as a function of 6 (at z = 6.5, upper panel), 
as a function of mean IGM neutral fraction ihi ( a * = 10', central panel), and as a function of redshift (at 8 = 10', lower panel). Right: 
Corresponding examples of skewness functions. 



tion there is suppressed, slowing down the subsequent ioni- 
sation rate relative to more neutral regions. Comparison of 
the fiducial model having (T min ,T ion ) = (10 4 K, 10 5 K) with 
the case of no feedback (T min ,T ion ) = (10 4 K, 10 4 K) illus- 
trates the suppression of both the amplitude of the auto- 
correlation and skewness functions in the presence of radia- 
tive feedback. The redshift where the minimum values of 
auto-correlation and skewness are reached is raised by the 
inclusion of feedback (given a common redshift for the com- 
pletion of reionisation) . The signature of radiative feedback 
is seen in the evolution of auto-correlation and skewness 
amplitude with neutral fraction. While there is substantial 
difference in the predicted amplitudes at neutral fractions 
below 0.5, the two cases are very similar before reionisation 
reaches a significant level. Thus the signature of radiative 
feedback is to suppress the amplitude of fluctuations only 
late in the reionisation epoch. 

In a model where the minimum mass for galaxy for- 



mation is 10 5 K in all regions (T min = T ion = 10 5 K), the 
amplitudes of both the auto-correlation and skewness are 
substantially increased (relative to the fiducial model) at 
all redshifts following the onset of reionisation. Similarly, a 
model with T m i n = 250K and Tl on = 10 J K predicts ampli- 
tudes of both the auto-correlation and skewness that are 
substantially smaller than the fiducial model. Indeed the 
variation in amplitude of the auto-correlation at fixed red- 
shift reaches an order of magnitude over the astrophysically 
plausible range of minimum masses. This prediction of the 
variation of fluctuation amplitude with the mass of ionising 
sources is consistent with the numerical simulations of Lidz 
et al. (2006). The variation is due too the larger (smaller) 
bias of the ionising sources, which leads to a larger (smaller) 
variation in ionisation state between over-dense and under- 
dense regions. This larger (smaller) variation in turn leads to 
a larger (smaller) amplitudes for both the auto-correlation 
and skewness functions. 
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Figure 6. Contour plots showing the evolution of auto-correlation (upper panels) and skewness (central panels) functions with redshift 
and angle, assuming a clumping factor of C = 10. Two models for the minimum mass of galaxy formation are shown, (T m i n , Ti on ) = 
(10 4 K, 10 5 K) [left panels] and (T m ; n , Tj on ) = (250K, 10 5 K) [right panels]. The keys to the right of the panels show the numerical values 
assigned to the various contours. The lower panels show the evolution of average neutral fraction for these models. 



All models in Figure [5] show minima of the variance and 
skewness at redshifts of z ~ 8 — 9. The duration and depth of 
these minima are sensitive to the bias of the ionising sources. 
In particular, the ratio of the auto-correlation amplitude 
at high redshift to the auto-correlation amplitude at the 
minimum varies between a factor of 10 and a factor of 100 
across the models considered. Similar results are found for 
the skewness. The redshift of the minimum is also sensitive 
to the bias of the ionising sources, though the dependence 
is not monotonic. On the other hand there is a monotonic 
relation between the neutral fraction at the redshift of the 
minimum and the bias, with more massive ionising sources 
generating a minimum at a smaller ionised fraction. Models 
with low clumping and models with large source mass both 
increase the amplitude of the auto-correlation and skewness. 
This degeneracy between models may be broken by looking 
at the location of the minima in variance and skewness as 
functions of redshift and neutral fraction. Assuming a fixed 
epoch of reionisation, both large clumping factors and mas- 
sive ionising sources delay the minima in time. Conversely, 



a large clumping factor results in the minima occurring at 
a smaller average neutral fraction, while massive ionising 
sources lead to the minima occurring at a higher neutral 
fraction. However in practice the neutral fraction will be 
very hard to measure due to instrumental and foreground 
effects, which may limit the utility of this effect. 

In summary the evolution of the variance and skew- 
ness in 21cm intensity maps with redshift and mean neutral 
fraction will probe physical quantities like the clumping of 
the IGM and the mass of ionising sources. In particular, 
reionisation scenarios with smaller values of clumping fac- 
tor, or more massive (and therefore more biased) ionising 
sources result in 21cm intensity maps with larger values of 
variance and skewness. Indeed, the variation in amplitude 
of the auto-correlation at fixed redshift reaches an order of 
magnitude over the astrophysically plausible range of min- 
imum masses. In addition, the relative locations of minima 
in variance and skewness when plotted as functions of red- 
shift and neutral fraction can be used to distinguish between 
amplitudes that are increased due to large values of ionising 
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source mass, or small values of clumping factor. While more 
massive sources increase the amplitude of fluctuations at all 
epochs, the signature of radiative feedback is the suppres- 
sion of fluctuations at epochs following, but not prior to the 
substantial reionisation of the IGM. We therefore conclude 
that the detection of fluctuations in 21cm intensity should 
unambiguously determine the mass of ionising sources. 



4 MEASUREMENT OF SKEWNESS IN 21CM 
MAPS 

In the previous sections we have demonstrated that galaxy 
bias should lead to skewness in the intensity distribution 
derived from redshifted 21cm maps of the reionisation epoch. 
In this section we discuss the prospects for measuring that 
skewness using first generation instruments. 

We first discuss the response of a phased array to the 
brightness temperature contrast of the IGM. Assuming that 
calibration can be performed ideally, and that foreground 
subtraction is perfect, the root-mean-square fluctuations in 
brightness temperature are given by the radiometer equation 



AT : 



(8) 



where A is the wavelength, T sys is the system temperature, 
Atot the collecting area, Qb the effective solid angle of the 
synthesised beam in radians, tint is the integration time, Av 
is the size of the frequency bin, and e is a constant that de- 
scribes the overall efficiency of the telescope. We optimisti- 
cally adopt e = 1 in this paper. 

As a concrete example we consider the case of the 
Mileura Wide Field Array-Low Frequency Demonstrator 
(hereafter LFEQ), though the measurement described could 
also be made using other instruments such as LOFAR0. In 
units relevant for the LFD, and at v — 200MHz, we find 
(Wyithe, Loeb & Barnes 2005) 



AT 



7.5 



(ii) mK 



At. 



( Av 



V 1MHz 



-1/2 



1.97 

Cbeam / \ j4.LFD 
1/2 



V lOOhr I 



(9) 



Here Alfd is the collecting area of a phased array consisting 
of 500 tiles each with 16 cross-dipoles [the effective collect- 
ing area of an LFD tile with 4x4 cross-dipole array with 
1.07m spacing is ~ 17 - 19m 2 between 100 and 200MHz 
(B. Correy, private communication; Bowman, Morales & He- 
witt 2005)]. The system temperature will be dominated by 
the sky and in the frequency range of interest has a value 
T sys ~ 280[(1 + 2)/7.5] 2 3 K. Av is the frequency range over 
which the signal is smoothed and #bcam is the size of the syn- 
thesised beam. The value of #beam can be regarded as the 
radius of a hypothetical top-hat beam, or as the variance of 
a hypothetical Gaussian beam. The corresponding values of 
the constant Cbeam are 1 and 1.97 respectively. 

Equation <(9j gives the uncertainty AT in the 21cm 
signal measured from a cylinder of IGM, with radius 6 



1 see http://www.haystack.mit.edu/ast/arrays/mwa/index.html 

2 http://www.lofar.org/ 



and depth Av. For comparison with our calculated auto- 
correlation and skewness, we modify this in order to es- 
timate the noise in the signal originating from a spheri- 
cal region of radius R, so that 6 = R/cIa and Avjv = 
2R{cdt/dz)- 1 {l + z)' 1 where v is the redshifted 21cm fre- 
quency. We find that consideration of spherical regions in- 
creases the noise relative to equation © by a factor of 
\/l578. 

In any one synthesised beam, the measurement noise 
will be far in excess of the 21cm signal. However as has been 
discussed previously (Bowman, Morales & Hewitt 2005), 
the determination of the power-spectrum may still be made 
since it can be constructed from a large number of indepen- 
dent measurements (each of signal-to-noise below unity). In 
practice one could determine the auto-correlation function 
(or power-spectrum) by dividing the data stream in half, and 
performing a cross-correlation (Max Tegmark, private com- 
munication). Hence defining T as the measurement signal, 
and AT as the thermal noise, we have 

<[(T - (T» + ATr][(T - (T» + ATi]) 
= ((T-(T)) 2 ) 

+((T - <T))ATi) + ((T - <T))AT 2 ) + (AT 1 AT 2 ), (10) 

where ATi and AT2 are drawn from Gaussian uncertainty 
distributions of width \/2AT. Here, with the exception of 
the first term (which is the variance), the cross-correlation 
results in terms that average to zero in the limit of large 
numbers of measurements. In a finite experiment, the sum 
of these terms is the measurement noise. 

To estimate the noise in the measurement, we have 
made many Monte-Carlo realisations of the observed dis- 
tribution. For each measurement we obtain N values of ob- 
served brightness temperature by adding telescope noise to 
the intrinsic T drawn at random from the theoretical distri- 
bution. The error for each measurement can then be calcu- 
lated from the noise terms in the cross-correlation, and the 
uncertainty in the measurement estimated as the variance 
among the values of error. For a field of angular size #fi e id and 
a Gaussian synthesised beam of radius 6, the number of in- 
dependent pointings is approximated by N = [#fi e id/(2#)] 2 , 
where the 2 in the denominator reduces the effect of over- 
lapping synthesised beams. We assume #fi c id = 10° for the 
LFD. 

Comparisons of signal and noise are plotted in the up- 
per panels of Figure [7] In both the left and right panels, 
the auto-correlation is plotted as a function of angle assum- 
ing C — 10, and at z = 6.5, z — 8 and z = 10.0 (thin 
lines). For comparison we plot the estimate of noise (thick 
lines), assuming integrations of 100 hours (left panel) and 
1000 hours (right panel). We find that the auto-correlation 
function can be detected at redshifts below the minimum in 
lOOhr. In 1000 hours, higher redshift fluctuations can also 
be measured. 

Our signal to noise estimate is approximate, and does 
not account for several important issues. For example, the 
LFD will have a very compact antenna arrangement. As a 
result the noise is not independent in all synthesised beams. 
This is because short baselines only give independent in- 
formation on length scales associated with their separation. 
Therefore with a compact arrangement, it is only at very 
large scales that the pixels become independent. Moreover, 
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we have only considered a single frequency slice, whereas 
fluctuations can be measured in many slices at slightly dif- 
ferent redshifts. This latter effect will tend to increase the 
sensitivity achievable relative to our estimate, while the for- 
mer will tend to reduce it. 

In analogy with the above procedure one could deter- 
mine the skewness function by dividing the data stream in 
three, and performing a cross-correlation. Hence 

<[(T - (T» + ATi][(T - (T» + AT 2 ][(T - (T» + AT 3 ]} 
= <(T-<T» 3 ) 

+((T - <T)) 2 AT!) + <(T - <T)) 2 AT 2 ) + ((T - <T)) 2 AT 3 } 

+((T - <T))ATi AT 2 ) + <(T - <T))ATiAT 3 ) 

+((T- <T))AT 2 Ar 3 ) + (ATiAr 2 AT 3 }, (If) 

where ATi, AT 2 and AT3 are drawn from Gaussian uncer- 
tainty distributions of width %/3AT. As before, with the 
exception of the first term (which is the skewness), the 
cross-correlation results in terms that average to zero in the 
limit of large numbers of measurements. In analogy with the 
above, we have made many Monte-Carlo measurements of 
the error from the noise terms in the cross-correlation, and 
estimated the uncertainty in the measurement of skewness 
form the variance among these Monte-Carlo errors. 

Comparisons of signal and noise for the skewness are 
plotted in the lower panels of Figure [7] In both the left and 
right panels, the skewness is plotted as a function of angle 
assuming C = 10, and at z = 6.5, 2 = 7.0, z — 8 and 
2 = 10.0 (thin lines). These curves correspond to the auto- 
correlation curves described above. For comparison we again 
plot the estimate of noise (thick lines) , assuming integrations 
of 100 hours (left panel) and 1000 hours (right panel). In 100 
hours, the skewness could only be (marginally) detected at 
lower redshifts, and at angles near 9 ~ 12'. However in 1000 
hours, the skewness can be detected over a range of angles 
at redshifts below the minimum. 

Estimates of the signal-to-noise (SN) for measurements 
of the auto-correlation and skewness were obtained by taking 
the ratio of the signal to the root-mean-square of the sum 
of the noise terms in equations (|10lll|l . The upper panels 
of Figure [S] show estimates of SN as a function of 9 for 
integrations of 100 hours and 1000 hours using the LFD. 
The thin lines show SN for the auto correlation function, 
while the thick lines show results for the skewness. In 1000 
hours and at low redshift, the skewness can be detected by 
the LFD at a signal-to-noise as large as SN ~ 10. The SN 
for detection of skewness is a factor of several to 10 below 
that achievable for the auto-correlation function. 

The above results suggest that the skewness in red- 
shifted 21cm maps should be just detectable with the LFD. 
However a hypothetical follow up to the LFD would com- 
prise ~ 10 times the collecting area, with baselines of ~ 10s 
to 100s of kilometers (we refer to this as the MWA-5000). 
The signal-to-noise achievable with the MWA-5000 is plot- 
ted as a function of 9 for integrations of 100 hours and 1000 
hours. In 100 hours, an MWA-5000 could detect skewness at 
high signal to noise outside the minimum redshift. In 1000 
hours an MWA-5000 could detect skewness at all redshifts, 
and map out its evolution in detail. Detection of skewness 
would therefore be an important scientific driver for a second 
generation array. 



5 CONCLUSION 

Using a simple semi-analytic model for density dependent 
reionisation, we have quantified the extent to which galaxy 
bias induces skewness in the intensity distributions of red- 
shifted 21cm maps. Our model demonstrates that skewness 
of the distribution is present on all scales, and during the 
entire reionisation era. In particular, skewness is present on 
scales where no single region is fully reionised. These results 
are consistent with the detailed numerical modeling previ- 
ously presented by Lidz et al. (2006) . The skewness (like the 
variance) decreases with increasing angular scale, but is not 
monotonic with redshift or neutral fraction, even in a sim- 
ple reionisation scenario. A feature of our density dependent 
reionisation model is that it predicts a minimum in both the 
variance and skewness at the time when the ionisation frac- 
tion is of order 50%. This minimum is a generic feature of 
models that include galaxy bias, and is due to the transition 
between the high redshift era when the 21cm fluctuations 
are generated by the density field, and the lower redshift era 
when fluctuations are dominated by biased reionisation of 
the IGM. 

The large scale skewness described in this paper is re- 
lated to, but is different from from the excess power (or 
bump) induced in the power-spectrum at small scales fol- 
lowing the appearance of bubbles (Furlanetto et al. 2004). 
We are unable to describe the statistics of small scale fluctu- 
ations which require full numerical modeling (see e.g. Lidz 
et al. 2006). Skewness appears in regions of low to moderate 
average ionisation, and is a result of galaxy bias arising from 
galaxy formation within large scale over-densities. We note 
that if the ionisation fraction was independent of density, or 
even if it was a linear function of over-density, then there 
would be no skewness. Indeed the relation between ionisa- 
tion fraction and density would be linear in the absence of 
galaxy bias (and of recombinations). However galaxy bias 
ensures a non-linear relation, and as a result, skewness of 
the 21cm intensity maps will be present throughout reioni- 
sation, and on all scales. 

We have discussed the evolution of the amplitudes of 
the auto-correlation and skewness of 21cm intensity maps, as 
well as their evolution with redshift and mean IGM neutral 
fraction. The evolution of the variance and skewness in 21cm 
intensity maps will probe quantities like the clumping of the 
IGM, the mass of ionising sources and the presence of radia- 
tive feedback. For example, we find that reionisation scenar- 
ios with smaller values of clumping factor, or more massive 
(and therefore more biased) ionising sources result in 21cm 
intensity maps with larger values of variance and skewness. 
Indeed the variation in amplitude of the auto-correlation at 
fixed redshift reaches an order of magnitude over the astro- 
physically plausible range of minimum masses. Moreover, 
the relative locations of minima in variance and skewness 
when plotted as functions of redshift and neutral fraction 
can distinguish between amplitudes that are increased due 
to large values of ionising source mass, or small values of 
clumping factor. In addition, while more massive sources in- 
crease the amplitude of fluctuations at all epochs, the signa- 
ture of radiative feedback is the suppression of fluctuations 
at epochs following, but not prior to the substantial reioni- 
sation of the IGM. We therefore conclude that the detection 
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of fluctuations in 21cm intensity should therefore unambigu- 
ously determine the mass of ionising sources. 

We have estimated the sensitivity of the Mileura Wide- 
field Array Low-Frequency Demonstrator to skewness in 
21cm intensity maps, and find that the predicted level should 
be detectable in long integrations (~ 1000 hours) at redshifts 
nearing the end of reionisation. However an instrument with 
10 times the collecting area of the LFD could detect skew- 
ness across a range of redshifts. Skewness in 21cm maps, and 
a minimum in the variance during the reionisation era are 
generic predictions of models of reionisation that account 
for biased galaxy formation. The presence of this minimum 
in variance, and of skewness, which will be detectable at a 
signal-to-noise a few to ten times smaller than the variance 
will provide confirmation of the cosmological nature of any 
observed angular fluctuations. 
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